home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / nrpas13.zip / COSFT.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  1KB  |  56 lines
  1. PROCEDURE cosft(VAR y: glyarray; n,isign: integer);
  2. (* Programs using routine COSFT must define the type
  3. TYPE
  4.    glyarray = ARRAY [1..n] OF real;
  5. where n is the dimension of the input data array. *)
  6. VAR
  7.    enf0,even,odd,sum,sume,sumo,y1,y2: real;
  8.    theta,wi,wr,wpi,wpr,wtemp: double;
  9.    jj,j,m,n2: integer;
  10. BEGIN
  11.    theta := 3.14159265358979/n;
  12.    wr := 1.0;
  13.    wi := 0.0;
  14.    wpr := -2.0*sqr(sin(0.5*theta));
  15.    wpi := sin(theta);
  16.    sum := y[1];
  17.    m := n DIV 2;
  18.    n2 := n+2;
  19.    FOR j := 2 TO (m+1) DO BEGIN
  20.       wtemp := wr;
  21.       wr := wr*wpr-wi*wpi+wr;
  22.       wi := wi*wpr+wtemp*wpi+wi;
  23.       y1 := 0.5*(y[j]+y[n2-j]);
  24.       y2 := (y[j]-y[n2-j]);
  25.       y[j] := y1-sngl(wi)*y2;
  26.       y[n2-j] := y1+sngl(wi)*y2;
  27.       sum := sum+sngl(wr)*y2
  28.    END;
  29.    realft(y,m,+1);
  30.    y[2] := sum;
  31.    FOR jj := 2 TO m DO BEGIN
  32.       j := 2*jj;
  33.       sum := sum+y[j];
  34.       y[j] := sum
  35.    END;
  36.    IF (isign = -1) THEN BEGIN
  37.       even := y[1];
  38.       odd := y[2];
  39.       FOR jj := 1 TO (m-1) DO BEGIN
  40.          j := 2*jj+1;
  41.          even := even+y[j];
  42.          odd := odd+y[j+1]
  43.       END;
  44.       enf0 := 2.0*(even-odd);
  45.       sumo := y[1]-enf0;
  46.       sume := (2.0*odd/n)-sumo;
  47.       y[1] := 0.5*enf0;
  48.       y[2] := y[2]-sume;
  49.       FOR jj := 1 TO (m-1) DO BEGIN
  50.          j := 2*jj+1;
  51.          y[j] := y[j]-sumo;
  52.          y[j+1] := y[j+1]-sume
  53.       END
  54.    END
  55. END;
  56.